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Abstract 

Following  the  work  of  Roberts  (Los  Alamos  Scientific  Laboratory  Report  No.  LA-299,  1945), 
we  investigate  the  effect  of  small  two-dimensional  perturbations  on  an  isolated,  planar  shock  front 
moving  steadily  through  an  inviscid  fluid  medium  with  an  arbitrary  equation  of  state  (EOS).  In 
the  context  of  an  initial-value  problem,  we  derive  explicit  analytical  expressions  for  the  linearized, 
time-dependent  Fourier  coefficients  associated  with  an  initial  corrugation  of  the  front.  The  tem¬ 
poral  evolution  of  these  coefficients  superficially  resembles  the  attenuated  “ringing”  of  a  damped 
harmonic  oscillator,  but  with  the  important  distinctions  that  the  frequency  of  oscillation  is  not 
constant,  and  that  the  damping  factor  is  not  simply  an  exponential  function  of  time  t.  It  is  shown 
that  at  least  two  three-parameter  families  of  stable  solutions  exist,  one  more  strongly  damped 
than  the  other.  In  both  cases,  we  find  that  the  envelope  of  oscillations  decays  asymptotically  as 
f~3/2,  with  shorter  wavelengths  dying  out  earlier  than  longer  ones.  For  a  particular  perturbed- 
shock  system,  the  strength  of  the  front  and  the  EOS  properties  of  the  material  through  which  it 
propagates  determine  the  applicable  family  of  solutions.  Theoretical  predictions  agree  well  with 
FAST2D  numerical  simulations  for  several  examples  derived  from  the  CALEOS  library. 
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I.  INTRODUCTION 


Shock  waves  are  ubiquitous  in  compressible  hydrodynamic  systems  ranging  from  inter¬ 
stellar  media  [1]  to  traffic  flow  [2],  Of  fundamental  importance  to  the  dynamics  of  shock 
waves  is  an  understanding  of  how  disturbances  to  otherwise  steady  shock  fronts  evolve  in 
time.  The  earliest  analysis  of  this  problem  is  likely  due  to  Roberts  [3],  who  considered  the 
effect  of  small  perturbations  on  isolated  planar  shocks  propagating  through  homogeneous 
fluids  from  the  point  of  view  of  an  initial-value  problem.  The  use  of  the  word  “isolated” 
here  implies  that  a  shock  is  very  far,  or  “de-coupled,”  from  its  driving  mechanism  ( e.g .,  a 
moving  piston).  Working  in  the  context  of  an  ideal  fluid  description,  Roberts  reduced  the 
linearized  system  of  governing  equations  to  an  integral  expression,  the  evaluation  of  which 
yields  the  time-dependent  amplitude  of  each  Fourier  component  of  the  initial  disturbance. 
His  work  demonstrated  that  shocks  are  almost  always  stable,  with  perturbations  decaying 
asymptotically  in  time  t  at  least  as  fast  as  f-1/2.  One  obvious  oversight  of  Roberts’  anal¬ 
ysis,  though,  is  that  the  presence  of  entropy  perturbations  behind  the  shock  is  neglected; 
moreover,  results  are  specialized  to  the  case  of  a  perfect  gas. 

In  this  paper,  we  generalize  Roberts’  calculation  and  derive  explicit  expressions  governing 
the  temporal  evolution  of  a  rippled  shock  front  in  a  fluid  medium  with  an  arbitrary  equation 
of  state  (EOS).  The  solution  methodology  is  based  on  a  perturbative  expansion,  with  the 
ratio  of  shock  ripple  amplitude  to  wavelength  serving  as  a  small  parameter.  Information 
about  the  EOS  of  a  specific  material  is  contained  in  three  dimensionless  parameters  that 
characterize  the  unperturbed  shock,  and  serve  as  necessary  inputs  for  the  first-order  theory. 
The  three  parameters,  which  in  this  study  are  derived  from  the  CALEOS  library  [4]  at  the 
U.S.  Naval  Research  Laboratory  (NRL),  are:  the  D’yakov  parameter  h  (a  quantity  that  is 
inversely  proportional  to  the  slope  of  the  shock  Hugoniot  [5]),  the  compression  ratio  r/,  and 
the  Mach  number  behind  the  shock,  Ad\.  Several  examples  are  considered  in  our  analysis 
that  underscore  a  somewhat  unexpected  result.  It  is  shown  that  at  least  two  families  of 
damped  oscillatory  solutions  exist  for  initial  disturbances  localized  at  the  shock  front.  The 
determination  of  which  family  applies  to  a  particular  perturbed-shock  system  depends  on  the 
sign  of  a  certain  dimensionless  quantity  A,  which  is  an  algebraic  function  of  the  three  shock 
parameters  h,  77,  and  M\  only.  The  properties  of  the  two  families  of  solutions  are  similar  in 
that  they  share  the  same  late-time  asymptotic  behavior  (oscillations  decay  in  time  as  t-3/2), 
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but  differ  in  the  degree  of  damping  that  the  oscillations  experience  initially;  see  Fig.  1.  All 
ideal  gases  belong  to  the  family  of  solutions  with  A  >  0.  For  moderately-strong  shocks  in 
materials  such  as  polystyrene,  aluminum,  and  deuterium-tritium  “ice,”  the  parameter  A  is 
negative  (according  to  the  CALEOS  database),  and  the  behavior  is  accurately  described  by 
the  more  strongly-damped  family  of  solutions.  At  sufficiently  large  driving  pressures,  though, 
the  sign  of  A  eventually  becomes  positive  —  the  threshold  for  this  continuous  crossover  being 
material  dependent. 

The  motivation  for  this  study  is  borne  out  of  an  interest  to  understand  better  the  implo¬ 
sion  of  inertial  confinement  fusion  (ICF)  targets,  which  contain  materials  whose  equations 
of  state  are  far  from  ideal  [6].  In  most  ICF  schemes,  it  is  required  to  compress  fuel  pellets 
containing  mixtures  of  deuterium  and  tritium  to  densities  as  much  as  1000  times  greater 
than  solid  values  [7] .  The  compression  of  the  fuel  results  from  multiple  shock  waves  launched 
into  the  pellet  by  high-intensity  radiation  striking  its  surface  and  ablating  away  the  outer¬ 
most  layers  of  material.  Successful  compression  requires  that  the  ablation  and  compression 
processes  occur  with  near-perfect  symmetry,  otherwise  high  fusion-reaction  yield  may  not  be 
achieved  [8].  Asymmetric  irradiation  and/or  rough  surface  finishes  can  lead  to  a  nonuniform 
ablation  process,  and  thus  to  the  generation  of  perturbed  shock  waves.  The  presence  of  dis¬ 
torted  fronts  in  ICF  pellets  is  significant  because  of  their  potential  for  seeding  hydrodynamic 
instabilities  (via  “interface  imprinting”  [9]  or  otherwise  [10]),  which  disturb  uniform  high- 
density  compression  and  reduce  gain.  A  better  understanding  of  the  dynamics  of  perturbed 
shocks  in  real  materials  could  provide  insight  into  how  to  suppress  these  instabilities,  thus 
improving  the  overall  uniformity  of  the  compression  process. 

Previous  work  by  Ishizaki  and  Nishihara  [11]  explored  the  subject  of  perturbed  shocks 
generated  by  nonuniform  ablation  surfaces  in  planar  ICF  targets.  The  authors  were  able 
to  solve  for  the  time  history  of  the  rippled  shock  and  ablation  fronts,  and  demonstrate  fair 
agreement  with  experiments  [12],  but  their  theoretical  model  relied  on  a  perfect-gas  EOS 
with  an  artificially  large  value  for  the  ratio  of  specific  heats.  In  the  present  study,  we  do  not 
attempt  to  model  the  complete  ICF-compression  scenario,  but  rather  focus  our  efforts  on 
assessing  the  influence  that  real  equations  of  state  have  on  the  dynamics  of  perturbed  shocks. 
Although  the  issue  of  how  such  perturbations  might  arise  physically  is  not  addressed  by  the 
theory,  we  expect  that  the  solution  contained  herein  will  contribute  to  a  more  complete 
understanding  of  the  shock-compression  process  in  ICF.  Moreover,  the  solution  may  have 
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direct  relevance  to  other  perturbed  shock  problems  in  gas  dynamics  as  well.  Fraley  [13],  for 
example,  has  shown  that  isolated  rippled  shocks  possess  many  of  the  same  general  properties 
as  those  generated  by  curved  pistons  [14,  15],  or  reflected  from  corrugated  walls  [16]. 

As  it  stands,  the  isolated  rippled-shock  solution  presented  in  this  paper  also  serves  as  a 
useful  benchmark  for  hydrodynamic-based  computer  codes  used  in  ICF  research.  Modern 
ICF  codes  are  usually  comprised  of  a  sequence  of  distinct  “modules”  for  simulating  vari¬ 
ous  physical  phenomena  such  as  hydrodynamic  motion,  thermal  conduction,  and  radiation 
transport,  as  well  as  accounting  for  realistic  EOS  information  (often  through  a  table  look-up 
procedure)  [17-21],  The  only  way  to  verify  the  fidelity  of  such  complicated  codes  is  through 
the  use  of  limiting  test  cases  that  are  amenable  to  analytical  solutions.  Presently,  however, 
very  few  analytically  tractable  problems  relevant  to  hydrodynamic  stability  calculations  are 
available  —  Rayleigh- Taylor  [22]  and  Richtmeyer-Meshkov  [23]  problems  being  the  only  other 
known  examples.  The  isolated  rippled-shock  problem  thus  provides  a  valuable  addition  to 
the  suite  of  available  “reality  checks”  for  testing  and  debugging  numerical  simulations.  This 
subject  has  been  discussed  previously  by  Munro  [24],  who  used  the  analytical  solution  of  a 
rippled  shock  problem  to  check  the  accuracy  of  the  code  LASNEX  [17],  but  found  rather 
poor  agreement.  In  contrast,  it  will  be  shown  here  that  the  FAST2D  code  [18,  21]  (with 
the  thermal  conduction  and  radiation  transport  modules  disabled)  performed  well  on  the 
isolated  rippled-shock  problem;  simulation  results  show  excellent  agreement  with  theoretical 
predictions  for  the  evolution  of  shock  ripple  amplitudes. 

This  paper  is  organized  as  follows.  In  Sec.  II,  we  formulate  the  initial-value  problem  for 
the  isolated  perturbed-shock  system,  and  derive  the  governing  first-order  equations,  as  well 
as  the  linearized  boundary  conditions  that  they  must  obey.  In  Sec.  Ill,  these  equations  are 
reduced  to  an  integral  expression  for  the  time-dependent  shock-ripple  amplitude.  (Up  to 
this  stage,  our  approach  closely  follows  that  of  Roberts,  but  with  the  important  distinctions 
that  the  present  analysis  accounts  for  the  existence  of  post-shock  entropy-vortex  waves,  and 
non-ideal-gas  equations  of  state.)  Section  IV  outlines  the  solution  of  the  integral  equation 
via  Laplace  transforms,  which  involves  one  of  two  different  inversion  procedures,  depending 
on  the  sign  of  the  quantity  A.  Comparisons  of  results  from  numerical  simulations  using 
the  FAST2D  code  in  conjunction  with  the  CALEOS  database  appear  in  Sec.  V.  Finally  in 
Sec.  VI,  the  conclusions  of  the  paper  are  given.  For  convenience,  Laplace  transforms  and 
other  mathematical  relations  that  were  useful  in  this  investigation  arc  listed  in  the  Appendix. 
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II.  FORMULATION  OF  THE  INITIAL- VALUE  PROBLEM 


The  first  step  in  our  analysis  is  to  define  the  “zeroth-order”  state  upon  which  perturba¬ 
tions  are  to  be  imposed.  We  choose  this  to  be  a  one-dimensional,  planar  step  shock  moving 
in  the  laboratory  frame  with  a  constant  speed  D.  The  shock  propagates  into  a  half-space 
filled  by  a  fluid  medium  at  rest  with  density  p0,  pressure  p0,  and  entropy  per  unit  mass  s0. 
Behind  the  unperturbed  shock,  the  density,  pressure,  and  specific  entropy  have  the  values 
Pi,  pi,  and  si,  respectively.  We  imagine  that  a  steadily-moving  driving  mechanism  (such 
as  a  piston  or  ablation  front)  supports  the  shock  from  behind,  but  that  its  influence  on  the 
shock  dynamics  is  negligible.  That  is,  we  assume  that  the  shock  is  far  enough  away  from 
the  driving  mechanism  that  the  transit  time  of  a  sound  wave  between  the  two  surfaces  is 
much  longer  than  the  time  of  interest  for  this  problem. 

The  unperturbed  state  (pi,pi)  is,  of  course,  not  arbitrary,  but  constrained  to  lie  on  the 
principal  Hugoniot  curve,  which  is  the  locus  of  all  compressed,  or  “downstream,”  states  that 
can  be  realized  behind  a  single  shock  front  given  the  initial,  or  “upstream,”  condition  (po,  po). 
The  shape  and  other  important  geometric  properties  of  this  curve  are  greatly  influenced 
by  the  EOS  of  the  substance  through  which  the  shock  propagates  [25].  In  the  case  of  a 
perfect  gas,  the  equation  for  the  Hugoniot  can  be  calculated  analytically,  with  the  result 
expressible  in  a  simple  closed  form.  For  non-ideal  materials,  the  situation  is  slightly  more 
involved  in  that  an  accurate  calculation  of  the  Hugoniot  usually  requires  the  use  of  complex 
EOS  models  to  supply  realistic  thermodynamic  data.  Examples  of  such  models  include 
CALEOS  [4],  SESAME  [26],  and  QEOS  [27],  all  of  which  characterize  a  variety  of  substances 
over  a  wide  range  of  pressure  and  density  states  by  combining  theoretical,  empirical,  and 
phenomenological  descriptions  [6].  A  few  Hugoniots  derived  from  the  CALEOS  database  are 
presented  in  Fig.  2  for  materials  of  interest  to  ICF  research.  For  comparison,  these  figures 
also  show  the  corresponding  results  for  the  more  widely-used  SESAME  library.  As  Fig.  2 
demonstrates,  the  theoretical  predictions  of  different  EOS  models  for  the  same  material  are 
not  always  in  agreement,  but  we  shall  not  attempt  to  address  these  discrepancies  here  since 
such  a  discussion  lies  beyond  the  scope  of  the  present  investigation.  In  the  analysis  that 
follows,  we  shall  employ  the  CALEOS  model  exclusively  to  characterize  the  properties  of 
shock  waves  in  materials  with  non-ideal  equations  of  state. 

A  quantity  that  plays  an  important  role  in  the  study  of  perturbed  shocks  is  the  slope  of 
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the  Hugoniot  curve  in  the  plane  of  density  p  versus  pressure  p.  Note  that  this  slope  —  which 
is  written  as  (dp/ dp)  h  and  evaluated  at  the  downstream  state  —  is  distinct  from  the  post¬ 
shock  isentropic  derivative,  (dp/dp)s  =  c2.  Here,  the  subscript  s  implies  constant  entropy, 
and  c  is  the  speed  of  sound  in  the  compressed  fluid.  Although  these  two  derivatives  possess 
a  second-order  tangency  at  the  initial  state,  and  are  approximately  equal  for  weak  shock 
waves  [29],  they  can  differ  appreciably  even  for  moderately  strong  shocks.  A  useful  way 
of  expressing  the  (inverse)  slope  of  the  Hugoniot  is  in  terms  of  the  dimensionless  “D’yakov 
parameter”  [5] 

k=-4 (t 

Pi  \dp 

where  j  =  p^D  is  the  mass  flux  density  across  the  shock.  According  to  a  linearized  normal¬ 
mode  analysis  [5,  30,  31],  stable  shocks  obeys  the  condition 

-1  <  h  <  1  +  2Mi ,  (2) 

where  Mi  =  (D  —  U)/c  is  a  downstream  Mach  number  satisfying  0  <  Mi  <  1,  and  U 
is  the  mass  velocity  of  the  shocked  fluid  in  the  laboratory  reference  frame.  Planar  shocks 
in  materials  with  values  of  h  that  lie  outside  of  this  range  experience  exponential  growth 
of  perturbations  [25,  32],  and  are  thus  unstable.  In  Eq.  (2),  the  lower  limit  is  known  to 
correspond  to  the  breakup  of  a  shock  into  two  waves  traveling  in  the  same  direction  [33], 
and  has  long  been  observed  experimentally  in  substances  undergoing  phase  transformations, 
or  yielding  plastically  at  the  clastic  limit  [34],  Satisfaction  of  the  condition  1  +  2Mi  <  h,  on 
the  other  hand,  has  been  shown  by  Gardner  [35]  to  correspond  to  the  splitting  of  a  shock 
into  two  counter-propagating  waves,  although  this  phenomena  apparently  lacks  experimental 
confirmation.  For  the  examples  considered  in  this  study,  all  values  of  h  lie  well  within  the 
stability  limits  specified  by  Eq.  (2). 

Let  us  proceed  with  our  perturbative  analysis  of  the  isolated  rippled-shock  problem.  In 
anticipation  of  the  linearization  procedure  to  follow,  it  is  convenient  at  this  stage  to  transform 
our  frame  of  reference  to  one  in  which  the  unperturbed  compressed  fluid  behind  the  shock 
front  is  stationary.  In  such  a  reference  frame,  the  unperturbed  planar  shock  moves  with  speed 
D—U ,  the  upstream  fluid  flows  at  speed  U,  and  the  normal  and  transverse  components  of  the 
unperturbed  velocity  field  behind  the  shock  front  vanish.  This  arrangement  is  advantageous 
because  it  reduces  the  complexity  of  the  first-order  algebraic  expressions  significantly.  By 
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applying  the  law  of  mass  conservation,  one  can  easily  show  that 


D  —  U  =  D/rj ,  (3) 

where  r)  =  pi/po  and  D  >  U.  Additionally,  conservation  of  momentum  requires 

Pi  =  Po  +  Pi  D2  (r)  -  l)/rj2  .  (4) 

Equations  (3)  and  (4)  are  two  of  the  three  well-known  conservation  laws  for  planar  shocks 
known  as  the  Rankine-Hugoniot  relations  [29] .  The  missing  equation,  which  expresses  con¬ 
servation  of  energy,  does  not  play  a  direct  role  in  the  present  analysis. 

Next,  we  wish  to  derive  an  expression  for  the  position  of  the  propagating  shock  front. 
To  do  so,  we  let  x  be  the  coordinate  normal  to  the  undisturbed  planar  front,  and  y  the 
coordinate  along  it,  as  shown  in  Fig.  3.  Assuming  that  the  shock  moves  in  the  negative 
^-direction,  its  unperturbed  position  as  a  function  of  time  is  given  by  x  +  (D  —  U)t  =  0.  We 
now  introduce  first-order  perturbations  into  the  shape  of  the  front,  and  in  the  hydrodynamic 
quantities  behind  it.  In  general,  such  a  disturbance  will  give  rise  to  two  types  of  waves  in  the 
downstream  flow:  entropy- vortex  and  sound  waves  [25] .  Since  these  waves  carry  energy  away 
from  the  front  whose  distortion  is  the  ultimate  source  of  the  perturbation,  we  may  expect 
that  it  will  eventually  regain  its  planarity,  even  in  the  absence  of  transport  coefficients  such 
as  viscosity  and  thermal  conduction  [36] .  For  a  single-mode  perturbation  with  wavenumber 
k  and  amplitude  5x(t),  the  position  of  the  perturbed  shock  is 

xs(y,t)  =  -(D  -  U)t  +  6x(t)elky,  (5) 

where  i  =  and  the  real  part  of  the  right  side  of  this  equation  is  implied.  Note  that 

to  account  for  an  arbitrary  multi-mode  perturbation  to  the  shock  front,  the  8x  in  Eq.  (5) 
should  be  replaced  with  the  indexed  coefficient  Sxk,  and  a  summation  over  all  /e-values  in  the 
spectrum  performed.  (A  similar  Fourier  decomposition  would  apply  to  the  hydrodynamic 
quantities  to  be  introduced  shortly.)  Owing  to  the  linear  independence  of  the  functions 
exp  (iky),  though,  the  result  of  this  procedure  is  merely  a  single  set  of  relations  that  must 
be  satisfied  for  every  value  of  k,  and  is  identical  to  that  obtained  using  the  single-mode 
formulation  above.  For  this  reason,  we  shall  avoid  the  additional  notational  complexity 
associated  with  a  formal  superposition  of  Fourier  modes,  and  limit  our  analysis  to  a  single 
h-value  only.  The  final  expressions  that  we  derive  can  then  be  easily  generalized  to  treat 
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multi-mode  perturbation  as  the  need  arises.  We  should  also  point  out  that  the  amplitude 
8x(t)  in  Eq.  (5)  is  considered  to  be  a  “small”  quantity  —  a  statement  that  will  be  made 
more  precise  momentarily.  Furthermore,  note  that  at  this  stage  the  functional  dependence 
of  8x(t )  is  unconstrained;  we  have  not  assumed,  for  instance,  that  it  has  the  “normal-mode” 
form  exp(— icut),  where  the  frequency  u  would  be  related  to  k  through  a  dispersion  relation. 

Since  the  zeroth-order  flow  is  stationary  in  the  shocked-gas  reference  frame,  we  can  write 


the  post-shock  normal  and  transverse  components  of  the  flow  velocity  as 

ux(x,y,t )  =  5ux(x,t)  elky  ,  (6) 

Uy(x,y,t )  =  5uy(x,  t)  elky  ,  (7) 

respectively.  Similarly,  the  pressure,  density,  and  specific  entropy  behind  the  shock  are 

p(x,  y,t)  =  Pi  +  5p(x,t)  elky  ,  (8) 

P(x,y,t )  =  pi  +  Sp(x,  t)  elky  ,  (9) 

s(x,y )  =  si  +  Ss(x)  elky  .  (10) 


In  Eqs.  (8)-(10),  a  quantity  prefixed  by  a  5  denotes  the  amplitude  of  a  hydrodynamic  per¬ 
turbation  whose  magnitude  is  assumed  to  be  much  smaller  than  its  zeroth-order  counterpart 
(e.g.,  8p  «  pi).  Note  that  in  the  shocked-gas  reference  frame,  the  entropy  perturbation  5s 
is  independent  of  time.  This  is  a  consequence  of  the  fact  that  entropy-vortex  waves  behind  a 
shock  propagate  with  the  zeroth-order,  downstream  fluid  velocity,  which  has  zero  magnitude 
here. 

Let  us  determine  the  first-order  equations  describing  hydrodynamic  motion  behind  the 
perturbed  shock  front.  The  governing  expressions  are  the  Euler  equations  for  isentropic  flow: 


dp 

dt 


+  V  ■  (pu)  =  0, 


chi  , 

p ~dt  +p('u'  V^u 


=  -vp, 


ds 

dt 


u-  Vs  =  0 . 


in) 

(12) 

(13) 


Here,  the  fluid  velocity  u  is  given  by  uxx  +  uyy,  where  x  and  y  are  unit  vectors  in  the  x-  and 
y-directions,  respectively.  Substituting  Eqs.  (6)-(9)  into  Eq.  (11),  neglecting  terms  involving 
products  of  perturbation  quantities,  and  canceling  common  exponential  factors,  we  arrive 


at  the  linearized  continuity  equation: 


d_ 

dt 


bp  +  pi 


d_ 

dx 


bux  +  ik  buv  1=0. 


(14) 


Similarly,  the  x-  and  ^-components  of  the  momentum  equation  reduce  to 

9  tr  ,  9  Sp  _  n 

—  OUx  +  77 - ~  0  , 

Ot  OX  Pi 

0  bp 

—  OUy  +IK  -  =  0  . 

Ot  pi 


(15) 

(16) 


In  deriving  Eq.  (14),  we  have  used  the  linearized  form  of  Eq.  (13),  i.e.,  d(Ss)/dt  =  0,  and 
the  thermodynamic  relation  5p  =  bp/c 2  +  ( dp/ds)p5s .  The  partial  derivative  with  respect 
to  time  of  this  latter  equation  allows  us  to  write 


0 


0 


mSp~c  mSp  =  0- 


(17) 


Equations  (14)-(17)  constitute  the  linearized  system  of  perturbed  fluid  equations  governing 
the  dynamics  of  a  rippled,  two-dimensional  shock  wave.  Before  they  can  be  solved,  though, 
these  equations  must  be  supplemented  by  boundary  conditions,  which  we  now  discuss. 

The  boundary  conditions  for  this  problem  are  derived  by  applying  the  principles  of  mass 
and  momentum  conservation  across  the  shock  front.  To  apply  these  principles,  we  first  need 
to  find  expressions  for  the  normal  and  tangential  unit  vectors  on  the  surface  of  the  rippled 
shock  (see  Fig.  3).  In  the  approximation  that  is  linear  in  the  perturbation  amplitudes,  these 
unit  vectors  are  given  by 


N  =  x  —  ik5xelky  y  , 
T  =  ikSx  elkyx  +  y  . 


(18) 

(19) 


Note  that  in  writing  Eqs.  (18)  and  (19),  we  have  implicitly  assumed  that  kSx  «  1  - 
a  statement  that  gives  precise  meaning  to  the  term  “small”  perturbation  in  the  present 
context.  With  us  =  xdxs/dt  denoting  the  perturbed  shock  velocity,  conservation  of  mass 
requires  that  p(us  —  u)  •  N  be  equal  on  both  sides  of  the  front.  This  leads  to  the  first-order 
equation 


5ux 


q  —  Id  D  bp 

'  0 «// 

77  dt  77  pi 


(20) 


The  zeroth-order  contribution  to  the  mass  conservation  equation  is  given  by  Eq.  (3). 
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We  now  consider  the  implications  of  momentum  conservation  across  the  rippled  shock. 
This  principle  states  that  the  vector  p(us  —  u)(us  —  u)  •  N  +  p  N  is  the  same  on  both  sides 
of  the  front,  which  to  zeroth-order  yields  Eq.  (4).  Using  the  relation 


c  |  dp\  D2 

Sp={TP)Ils^-Wh 


ip, 


(21) 


one  can  show  after  performing  some  algebra  that  the  x-component  of  the  first-order 
momentum-conservation  equation  can  be  written  as 


2D  pi 


(22) 


It  should  be  noted  that  the  relation  between  the  first-order  pressure  and  density  amplitudes 
in  Eq.  (21)  holds  immediately  behind  the  shock  front  only,  and  is  merely  a  consequence 
of  expanding  the  Hugoniot  equation  in  a  first-order  Taylor  series.  Equation  (22)  can  be 
combined  with  Eq.  (20)  to  yield 


Sux  = 
Sp  = 


1  f  1  —  h\  d 


Sx 


rj  \  1  +  h )  dt 
2D  pi  f  rj  —  1  \  d 


Sx . 


(23) 

(24) 


1  +  h  V  rf  J  dt 

To  complete  our  set  of  boundary  conditions,  we  also  need  an  expression  for  the  transverse 
amplitude  Suy.  This  is  obtained  by  enforcing  continuity  of  the  transverse  component  of 
velocity  across  the  shock  front.  The  result  is 


Suv  =  ikD  — - Sx  . 

V 


(25) 


The  presence  of  the  factor  %  in  Eq.  (25)  implies  that  the  transverse  velocity  perturbation  is 
90°  out-of-phase  with  the  ripple  on  the  surface  of  the  shock. 

In  order  to  simplify  the  statement  of  our  problem,  it  is  useful  at  this  stage  to  introduce 
the  normalized  space  and  time  variables  £  =  kx  and  r  =  Dkt/r],  along  with  the  following 
definitions: 

r)  —  1 


9(j)  =  k 

rj 

H(,r)  = 

D2  pi 

Tj 

Vm£>t)  =  p8ux, 


Sx , 


(26) 

(27) 

(28) 
(29) 
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Employing  Eqs.  (26)-(29),  the  linearized  mass  and  momentum  equations  [Eqs.  (14)-(16)]  can 
be  cast  in  dimensionless  form  as 


i'l/jy  + 


d^x  d<p_ 
dr  d£ 


di>y 

lh 


+  i(fi  =  0, 


(30) 

(31) 

(32) 


where  we  have  used  Eq.  (17)  to  eliminate  Sp  from  Eq.  (14).  These  equations  are  subject  to 
the  boundary  conditions 


t. 


h 


-T,  T  = 


1  +  h 

^y(-T,r)  =  ig  g(r) , 
2 

0(-r,r)  = 


1  +  h 


9\r) , 
9(r) , 


(33) 

(34) 

(35) 


at  the  shock  front,  whose  position  in  normalized  variables  is  given  by  £  =  — r.  Note  that  in 
Eqs.  (33)  and  (35),  the  prime  on  g(r)  denotes  differentiation  with  respect  to  the  variable  r. 

The  system  of  first-order,  partial  differential  equation  appearing  in  Eqs.  (30)-(32)  can 
be  reduced  to  a  single  second-order  equation  for  the  dimensionless  pressure  amplitude  <f>  by 
taking  the  r-derivative  of  Eq.  (30),  subtracting  from  it  the  ^-derivative  of  Eq.  (31),  and 
employing  Eq.  (32).  Following  this  procedure,  we  End 

d2(j) 


+  M[ 


2d V 


+  0  —  0  , 


(36) 


d £2  dr 2 

which  we  recognize  as  a  homogeneous,  two-dimensional  wave  equation  that  has  been  Fourier 
transformed  in  one  spatial  variable.  In  passing,  we  note  that  Eq.  (36)  is  equivalent  to 
the  Klein-Gordon  equation  for  describing  the  dynamics  of  a  “scalar”  meson  in  quantum 
mechanics  [37];  an  equation  with  the  same  form  as  Eq.  (36)  also  governs  the  behavior  of 
a  flexible  string  with  additional  stiffness  forces  provided  by  its  surrounding  medium  ( e.g ., 
a  string  embedded  in  a  sheet  or  rubber)  [38].  This  equation  is  subject  to  Cauchy-type 
(Dirichlet  and  Neumann)  conditions  at  the  initial  instant,  and  at  the  surface  of  the  shock. 
Since  Eq.  (36)  is  hyperbolic,  this  a  well-posed  problem  possessing  a  unique,  stable  solution 
[39].  The  initial  conditions  that  (j)  must  obey  are 


o(C0)  =  MZh 


!«•">  - 


-M[ 


-2 


#»K.o)  +  ^K,o) 


(37) 

(38) 
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where  the  latter  expression  results  from  Eq.  (30).  At  the  shock  front,  0  must  satisfy  Eq.  (35) 
at  all  times,  as  well  as  a  Neumann-type  boundary  condition,  which  can  be  derived  by 
subtracting  Eq.  (30)  from  Eq.  (31),  and  using  the  total  r-derivative  of  Eq.  (33),  where 


d  d  d 
dr  d  t  <9£ 


when  £  =  —  r. 


The  result  is 


h 


1  +  h 


g"{r)  +  vg(T)  = 


<90  2  <90 

<9£  1  <9r 


(39) 


Z=-T 


The  remainder  of  this  paper  is  devoted  to  finding  an  explicit  solution  to  Eq.  (36),  subject  to 
the  conditions  specified  in  Eqs.  (33)-(35)  and  Eqs.  (37)-(39).  We  should  point  out,  though, 
that  the  result  of  our  analysis  does  not  yield  the  normalized  pressure  amplitude  0  directly; 
instead  we  choose  to  solve  for  the  time-dependent  ripple  amplitude  g(r),  since  that  quantity 
is  more  easily  compared  with  results  from  numerical  simulations.  Once  g(r)  is  determined, 
the  function  0  can  in  principle  be  found  by  recourse  to  the  governing  system  of  equations 
(although  that  step  is  not  performed  in  the  present  study). 

It  is  somewhat  instructive  to  compare  Eq.  (39)  with  the  equation  governing  damped 
harmonic  motion  [40]: 


mr"(t )  +  nr(t)  =  —vr'(t) . 


(40) 


Here,  the  quantities  m,  k  and  v  are  constants,  and  the  function  r(t)  represents  some  time- 
dependent  amplitude  (e.g.,  the  position  of  a  mass  connected  to  a  spring  on  a  frictional 
surface).  From  this  comparison,  we  see  that  Eq.  (39)  may  be  viewed  as  a  harmonic  oscillator 
equation  with  a  complicated  damping  term  —  i.e.,  one  that  is  not  simply  proportional  to 
the  velocity  of  motion.  As  we  shall  see,  the  evolution  of  the  ripple  amplitude  g(r)  in  Eq.  (39) 
does  superficially  resemble  the  attenuated  “ringing”  of  a  damped  harmonic  oscillator,  but 
with  the  important  distinctions  that  the  frequency  of  oscillation  is  not  constant,  and  that 
the  damping  factor  is  not  simply  an  exponential  function  of  time. 


III.  AN  INTEGRAL  EXPRESSION  FOR  THE  RIPPLE  AMPLITUDE 

In  this  section,  we  outline  our  strategy  for  solving  Eq.  (36),  which  is  a  linear,  hyperbolic, 
partial  differential  equation.  Following  Roberts  [3],  we  begin  our  analysis  by  noting  that 
this  equation  possesses  the  characteristics 

£-£a  =  ±(ra-r)/M1.  (41) 
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Here,  £a  and  ra  are  constants,  and  refer  to  an  arbitrary  point  in  the  plane  define  by  the 
variables  £  and  r.  If  we  choose  this  point  to  correspond  to  the  position  of  the  shock  at  the 
time  Ta  so  that  £a  =  — ra,  we  see  that  only  one  of  the  characteristics  comes  from  the  shocked 
material  —  namely  the  one  preceded  by  a  “+”  sign  in  Eq.  (41).  Thus,  we  conclude  that 
conditions  at  the  front  are  influenced  only  by  what  occurs  in  the  region  G  defined  by 

t  >  0 , 

£  +  r  >  0, 

£  +  Ta  <  {ra  —  t)/Mi  . 

This  region  is  bounded  by  the  contour  C,  which  is  comprised  of  three  line  segments,  as  shown 
in  Fig.  4.  The  Erst  extends  along  the  £-axis  from  0  to  ra(M f1  —  1),  the  second  along  the 
characteristic  £  +  ra  =  (ra  —  r)/Mi,  and  the  third  along  the  shock- wave  path  £  +  r  =  0.  For 
convenience,  we  shall  refer  to  these  three  line  segments  as  Ci,  Cn,  and  Cm,  respectively. 
Note  that  in  the  dimensionless  units  of  this  problem,  the  speed  of  sound  is  given  by  l/M\. 

We  now  seek  a  solution  to  Eq.  (36)  within  the  region  G.  Since  Cauchy  boundary  are 
to  be  imposed  on  only  two  of  the  three  bounding  line  segments  —  he.,  along  Ci  and  Cm, 
but  not  along  the  characteristic  Cn  —  this  region  is  classified  as  an  “open  surface,”  and  so 
we  may  expect  to  find  a  unique  and  stable  solution  there  [39] .  We  proceed  by  employing  a 
method  due  to  Riemann  for  solving  hyperbolic  partial-differential  equations  (analogous  to 
the  theory  of  Green  functions  for  elliptic  operators)  [38].  This  method  relies  on  choosing 
an  appropriate  “Riemann  function”  /(£,  r  |  £',  r'),  which  in  our  case,  should  be  a  particular 
solution  of  Eq.  (36).  The  function  /  must  also  obey  the  reciprocity  relation:  /(£,  t  |  £',  t')  = 
|£,  r).  Here,  the  primed  coordinates  refer  to  an  arbitrary  point  in  the  space-time 
diagram  of  Fig.  4.  Such  a  function  can  be  easily  found  by  assuming  that  it  depends  only  on 
the  quantity 


Substitution  of  f(R)  into  Eq.  (36)  leads  to  a  Bessel  equation,  with  the  solution 

m  =  Jo  W'  -  t> w  -  (F  -  £)2  • 

The  symbol  Jq  in  this  expression  denotes  a  Bessel  function  of  order  zero.  Note  that  the 
Neumann  function,  which  is  the  second  solution  of  Bessel’s  equation,  is  disqualified  for  use 
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in  the  present  context  since  it  is  singular  at  R  —  0.  For  the  primed  coordinates  in  f(R),  we 
choose  f  =  —  ra  and  r'  =  ra. 

Let  us  demonstrate  how  the  Riemann  function  f(R)  can  be  used  to  to  reduce  Eq.  (36)  to 
an  integral  equation  for  the  ripple  amplitude  g(r).  The  first  step  is  to  define  the  operator 


Q  = 


d2  d2 


d& 


+  Mf 


dr 2 


and  introduce  the  vector 


V 


df 

df 


'D 


f  +  Mf 


where  f  and  t  are  unit  vectors  in  the  f  and  r  directions,  respectively.  An  important  property 
of  V  is  that  it  is  solenoidal:  V  •  V  =  f'Qf  —  <f>Q  f  =  0  .  Next,  we  compute  the  line  integral  of 
this  vector  around  the  contour  C  by  employing  the  two-dimensional  form  of  Gauss’  theorem: 


V- \d2x=  <b  V-ndl, 


IG 


(42) 


c 


where  h  is  an  outward-pointing  unit  vector  and  d£  is  an  infinitesimal  line  segment.  Then, 
using  Eq.  (42),  and  the  expressions  for  hdt  along  C7,  Cn,  Cm  listed  in  Table  I,  we  find 

/■(l— Ml)  Ta/Mi 


0  =  -Mf 


i o 


r  d(j>  df\ 

fTr~*Tr),  ^ 


Cl 


1 0 


d<p  df 

Dr  ^  dr 


I  rd<j)  df 


dr 


Cn 


Mf  f 


~  f 


df 

~dt, 


df 

df 


dr . 


Cm 


Values  for  /,  df  /dr,  and  df  /df  are  also  given  in  Table  I.  Substituting  these  expressions 
into  the  equation  above,  making  use  of  Eq.  (35)  and  Eqs.  (37)-(39),  and  noting  that 


d  d  Id 
~drr  ~  fdr  ~  ARdf 


when  f  +  Ta  =  (ra  -  t)/Mx  , 


leads  to 


2Mi 

1  +  h 


g\ra)  + 


1  -  h 
1  +  h 


g"{r)  +  r)g(r) 


J o  cx(Ta 


7-)]  dr  =  x(ja) , 


where 


a 


1  -  Mf 
Mf  ’ 


(43) 


(44) 
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and 


In  certain  limiting  cases,  the  complexity  of  this  last  equation  can  be  reduced  considerably. 
For  example,  if  we  consider  an  initial  perturbation  consisting  of  a  sinusoidal  deformation  to 
a  planar  shock  front  (with  unperturbed  hydrodynamic  quantities  behind  it)  we  have 

x(t*)  =  \  ^(°)  Maro) »  (46) 

which  is  valid  for  ra  >  0.  This  equation  results  from  setting  the  quantities  0o(O,  0), 

and  ^(£,0)  —  but  not  c?t/>x(£,  0)/<9£  —  equal  to  zero  for  £  >  0  in  Eq.  (45).  Note  that  for 
this  initial  condition,  the  function  y  is  discontinuous  at  ra  =  0,  where  it  assumes  the  value 
2Mig'(0)/(l  +  h). 

In  order  to  cast  Eq.  (43)  in  a  more  useful  form,  we  manipulate  it  according  to  the  following 
procedure.  First,  we  change  the  variable  of  integration  from  r  to  6.  Next,  we  substitute 
the  variable  r  for  ra.  Finally,  we  perform  an  integration  by  parts,  and  then  integrate  the 
entire  resulting  expression  with  respect  to  r.  The  result  of  all  these  operations  is  a  Volterra 
equation  of  the  second  kind  [41]  for  the  ripple  amplitude: 


For  the  case  of  a  planar  shock  front  initially  deformed  into  a  sinusoidal  shape,  we  have 

F(r)  =  y — ~  h  [2M,  +  (1  -  h)  J0  (a  r)]  .  (50) 

Note  that  F( 0)  =  g(0).  In  the  next  section,  we  show  how  the  integral  equation  in  Eq.  (47) 
may  be  solved  in  this  special  case  using  the  method  of  Laplace  transforms. 
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IV.  THE  SOLUTION 


The  difficulty  associated  with  an  integral  expression  like  the  Volterra  equation  in  Eq.  (47) 
is,  of  course,  the  appearance  both  inside  and  outside  of  an  integral  sign  of  the  unknown 
function  that  we  seek  [in  our  case  the  ripple  amplitude  g{r)\.  This  fact  greatly  complicates 
a  straightforward  attempt  at  finding  a  solution.  In  such  situations,  one  must  often  resort  to 
indirect  mathematical  approaches  such  as  integral  transform  methods  that  map  the  sought- 
after  function  from  one  space  (in  our  case  r)  to  another  in  which  a  solution  is  readily 
determined.  Often,  the  challenge  then  lies  in  the  inversion  of  the  mapping  procedure  to 
express  the  final  answer  in  terms  of  the  original  variable. 

One  such  approach  —  the  method  of  Laplace  transforms  [41]  —  is  particularly  well-suited 
to  the  present  class  of  problems.  A  fundamental  axiom  of  this  method  is  the  convolution 
theorem  [see  Eq.  (A9)  in  the  Appendix],  which  permits  us  to  convert  Eq.  (47)  from  an 
integral  equation  for  g(r)  into  an  algebraic  one  for  ^(s)  —  our  shorthand  notation  for  the 
Laplace  transform  of  g(r).  Throughout  our  discussion,  we  shall  use  the  subscript  L  to  denote 
the  Laplace  transform  of  a  function,  where  s  is  the  associated  transform  variable  (not  to  be 
confused  with  the  specific  entropy  introduced  in  Sec.  II);  see  the  Appendix.  Applying  the 
convolution  theorem  to  Eq.  (47)  and  solving  the  resulting  expression  for  ^(s),  we  find 

9M  =  Fl(8)/[1-Kl(s)]  . 


Limiting  our  attention  to  the  special  initial  condition  of  a  planar  shock  front  deformed  by  a 
sinusoidal  ripple,  explicit  expressions  for  Fl  and  Kjj  can  be  easily  derived  using  Eqs.  (A2)- 
(A5)  in  the  Appendix.  The  result  is  that  the  transform  of  the  ripple  amplitude  can  be 
written  as 


£l(s)  _  Vs2  +  a2  +  Ps 
g{ 0)  s  Vs2  +  a2  +  P  s2  +  T 

where  we  have  introduced  the  definitions 


(51) 


P  = 


1  -h 

2  Ad i  ’ 


r  =  (1  +  h)g 

2Mi 


(52) 

(53) 


Our  principal  task  in  this  section  is  to  invert  Eq.  (51),  and  thus  determine  the  time-dependent 
ripple  amplitude:  g(r)  =  £_1{c/,l(s)},  where  the  symbol  £_1  denotes  the  inverse  Laplace- 
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transform  operator.  Since  this  equation  involves  the  quotient  of  irrational  functions,  though, 
such  an  inversion  is  not  a  trivial  exercise,  and  requires  special  consideration. 

The  first  step  to  finding  the  inverse  Laplace  transform  of  <?l(s)/<?(0)  is  to  rationalize  the 
denominator  in  Eq.  (51),  which  yields 

gL(s)  _  ((32  -  l)s3  +  (f3T  -  a2)s  +  TVs2  +  a2 

0(0)  ~~  (/32  —  l)s4  +  (2/3T  —  a2)s2  +  T2  '  1  j 

To  proceed,  we  must  now  adopt  a  particular  solution  methodology.  One  obvious  approach 
is  to  seek  a  solution  of  Eq.  (54)  by  analyzing  its  poles,  and  then  computing  the  appropriate 
Bromwich  integral  using  well-established  methods  from  analytic  function  theory  [39].  This 
procedure  is  somewhat  involved,  however,  owing  to  the  square  root  term  in  Eq.  (54),  which 
necessitates  consideration  of  a  branch  cut  in  the  complex  plane.  A  slightly  more-attractive 
possibility  for  inverting  Eq.  (54)  is  to  factor  it  into  a  series  of  paired  multiplicative  expres¬ 
sions  (partial  fractions)  each  of  which  can  be  recognized  as  a  Laplace  transform  of  a  known 
function;  the  products  of  terms  in  this  series  can  then  be  inverted  through  use  of  the  convo¬ 
lution  theorem.  In  this  paper,  we  choose  to  pursue  the  latter  strategy.  In  so  doing,  we  must 
make  sure  that  the  result  has  no  imaginary  component,  since  g(r)  is  strictly  a  real  quantity. 
This  consideration  leads  to  two  ways  of  factoring  Eq.  (54),  which  in  turn  yields  two  families 
of  solutions. 

For  reasons  that  will  become  clear  shortly,  the  factoring  method  appropriate  for  a  given 
ripple-shock  system  is  determined  by  the  sign  of  the  quantity  A  defined  as 

A  =  a4-  A(3Ta2  +  4  T2  .  (55) 

Plotted  as  a  function  of  a2,  this  expression  forms  a  parabola,  as  shown  in  Fig.  5.  The  roots 
of  A  =  0  are  easily  shown  to  be 

4  =  2 T(f3  ±  y/^1)  • 

Note  that  the  inequality  A  <  0  is  consistent  with  the  value  of  a2  lying  between  and  a 2+. 
The  condition  A  >  0,  on  the  other  hand,  implies  that  a2  <  d2_  or  <  a2,  but  apparently 
only  the  former  inequality  is  physical.  In  the  derivation  that  follows,  we  shall  assume  that 
a2  never  exceeds  a2^  when  A  >  0.  Additionally,  we  require  that  /3  >  1  and  T  >  0  always 
-  conditions  that  appear  to  be  satisfied  for  most  “well-behaved”  equations  of  state.  Let  us 
now  discuss  the  derivation  of  the  solution  for  each  sign  of  A  separately. 
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A.  The  case  A  >  0 


If  the  sign  of  A  is  positive,  the  denominator  in  Eq.  (54)  can  be  written  as  (/ 3 2  —  1)  times 
the  product  [s2  +  (a  +  b)2][s2  +  (a  —  b)2],  where 

1 


(a  ±  b)  = 


2  ((32  -  1)  L 
and  a  and  b  are  real  constants  given  by 


2/3T  -a2T  \/«4  -  4/?r«2  +  4r2 


(56) 


'2/3r-a2 

O  —  \  I  ~rrr^ - 7T  + 


4(/?2  —  1)  2V/^TT’ 


b  =  -< 


i2pr  -  a2  _  r 
4(/52  -  1)  2a/ f32  -  1 


Since  A  >  0,  the  right  side  of  Eq.  (56)  is  a  real  quantity  whose  sign  is  positive  by  virtue  of 
the  inequality  a2  <  a2  and  the  fact  that  \2/3T  —  a2\  >  A1/2  for  (3 2  >  1.  Equation  (54)  can 
then  be  factored  as 


Cl 


9l(s)  =  1 

g(0)  s2  +  (a  +  b)2  [ s  +  /s2  +  a2 


+  C2s 


+ 


C3 


s2  +  (a  -  b)2  [s  +  a/s2  +  a2 


+  c4s 


(57) 


where  the  (real)  constants  ci,  C2,  c3,  and  c 4  are  given  by 

a2r 


Cl  = 


\J a4  -  4/ra2  +  4r2 


=  -c3, 


_  1  r  -  a2/2 

(1  _  2  ol4  -  A(3Ta2  +  4 r2 


=  1  —  C4  . 


Consulting  Eqs.  (A6)-(A8)  in  the  Appendix,  we  see  that  Eq.  (57)  is  in  the  form  of  products 
of  transforms  involving  trigonometric  and  Bessel  functions.  Using  the  convolution  theorem 
[Eq.  (A9)],  this  expression  can  be  inverted  to  give 

fi'(r)  ci  fT  .  u  Ji(oiz)  ,  / 

— — —  =  -  /  sin  (a  +  b)(r  —  z)\  - -  dz  _|_  C2  Cos(a  +  b)r 

g{ 0)  a  +  b  J o  az 

c3  fT  .  w  , , ,  Ji(az)  ,  , 

H - 7  /  sm[(a  —  o)(r  —  z)J  — - - az  +  c4cos(a  —  b)r  , 


a  —  b 


1 0 


where  the  symbol  Ji  denotes  a  Bessel  function  of  order  one.  In  the  case  that  a±b  >  a  (which 
holds  for  all  ideal  gases)  one  can  show  using  Eq.  (A10)  that  all  purely  oscillatory  terms  in 
the  above  expression  cancel.  If  a  +  b  and/or  a  —  b  are/is  less  than  a,  some  oscillatory  terms 
persist  and  stationary  perturbations  —  that  neither  grow  or  attenuate  in  time  —  result.  We 


18 


shall  not  discuss  this  unusual  phenomenon  further  here,  but  simply  remark  that  it  is  likely 
associated  with  the  D’yakov-Kontorovich  instability  of  shock  waves  [5,  42-45]. 

Assuming  a  ±b  >  a,  the  solution  for  A  >  0  can  be  written  after  some  manipulation  as 


ha4  -  A(3Ya2  +  4T2  J 0 


,  ■  u  u  ■  7  ,Ji[a{r  +  z)]  ,  ,  . 

(a  cos  az  sin  bz  —  b  sin  az  cos  bz) - - - - —  dz  ,  (58] 

a(r  +  z) 


where  we  have  used  the  fact  that  ...  dz  =  . . .  dz  —  Jt°°  . . .  dz.  Employing  Eq.  (A10) 

in  the  Appendix,  one  can  show  that  the  right  side  of  Eq.  (58)  has  the  correct  normalization 
(i.e.,  assumes  the  value  unity)  for  r  =  0.  Furthermore,  from  the  asymptotic  form  of  the  first- 
order  Bessel  function,  J\(q)  ~  yj2/{irq)  cos  [q— 3/(47t)]  for  q  — >  oo,  we  see  that  the  amplitude 
g(r)  undergoes  oscillations  that  die  out  as  t-3/2  late  in  time.  This  asymtotic  dependence, 
which  has  been  observed  previously  in  both  shock-tube  [16,  46]  and  laser-driven  ICF  [12] 
experiments,  is  apparently  a  general  property  of  perturbed  shock  fronts  that  extends  beyond 
the  “isolated”  variety  considered  in  this  paper  [13]. 


B.  The  case  A  <  0 


If  the  discriminant  A  is  negative,  Eq.  (54)  must  be  factored  differently  to  yield  a  real 
expression  for  g(r)/g( 0).  In  this  case,  we  find  that  Eq.  (57)  should  be  written  as 


9l(s) 

9(0) 


1  d\  -\-  (I2  s 

(• s  +  cr)2  +  a 2  _s  +  a/s2  +  a2 

i  r  d§  +  <^6 s 


T  d%  s  T  c?4 


(s  -  a)2  +  a2  Ls+  Vs2  +  a2 


+  d’j  s  +  ds  , 


(59) 


where  the  quantity  a  =  ib  is  real  and  positive.  The  constants  d\,  d^,  d:s ,  and  d^  in  Eq.  (59) 
are  given  by 
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Once  again,  Eq.  (59)  has  the  form  of  products  of  transformed  functions  that  can  be  easily 
recognized.  Using  the  convolution  theorem  and  Eqs.  (A6)-(A8)  in  the  Appendix,  we  find 
that  the  solution  for  A  <  0  is 


g(r) 

9(0) 


=  -e  <TT  cos  ar 


a 


of 


l(32  -  1  2(/32  -  1) 


4aa 


^\[W 


[  [  e-a(T-z) 

cos  air  —  z)  H —  sin  air  —  z) 

l  Jo 

L  a  J 

Ji(az) 

az 


dz 


+ 


a  . 

cos  az  H —  sm  az 
a 


J\  [a(r  +  z)] 
a(r  +  z) 


dz 


(60) 


In  arriving  at  this  expression,  we  have  made  use  of  Eqs.  (All)  and  (A12)  to  cancel  all 
non-evanescent  terms.  The  presence  of  decaying  exponential  functions  in  Eq.  (60)  tends  to 
enhance  the  initial  damping  of  shock-ripple  oscillations,  and  thereby  serves  to  distinguish 
this  family  of  solutions  from  that  in  Eq.  (58).  Note,  though,  that  for  both  families  the 
asymptotic  dependence  as  r  — »  oo  is  the  same,  namely  t-3//2  times  an  oscillatory  function 
of  r. 

The  same  asymptotic  behavior  holds  if  the  discriminant  A  vanishes  —  an  event  that  can 
occur  at  isolated  points  along  the  Hugoniot  curve  for  realistic  equations  of  state.  In  this 
case,  the  roots  (a±6)2  in  Eq.  (56)  are  identical  and  the  factoring  of  gjJ(s)/g( 0)  is  somewhat 
different  than  in  Eq.  (57)  or  Eq.  (59).  We  shall  not  provide  the  details  of  the  calculation  for 
A  =  0,  but  simply  quote  the  final  result,  which  is 
9(t) 

g(  0) 


=  V a2  —  b2  — 


cr 


^sin  V a2  —  b2  z  —  V a2  —  b2  z  cos  a Ja2  —  b2  z'j 


J\  [«(t  +  z)\ 
a(r  +  z) 


dz . 
(61) 

The  quantities  beneath  a  radical  sign  in  the  above  expression  are  positive  if  the  smaller  root 
of  Eq.  (55)  is  assumed  —  i.e.,  a2  =  a2  when  A  =  0.  It  should  be  emphasized  here  that 
the  ripple  amplitude  g{r)  undergoes  a  continuous  transition  from  one  family  of  solutions  to 
the  other  as  A  changes  sign.  That  is,  in  the  limit  that  A  — >  0,  both  Eq.  (58)  and  Eq.  (60) 
smoothly  approach  the  solution  appearing  in  Eq.  (61). 

We  should  also  remark  that  although  Eqs.  (58)  and  (60)  were  obtained  by  analyzing  a 
single  Fourier  mode  of  the  shock-front  perturbation,  the  same  solutions  apply  to  all  nor¬ 
malized  amplitudes  in  a  linearized  multi-mode  description.  Since  g(r)  and  g( 0)  are  both 
proportional  to  k,  the  wavenumber  enters  the  normalized  solution  g(r)/g( 0)  only  through 
the  independent  variable  r  =  Dkt/rj.  As  a  result,  we  see  that  shorter-wavelength  perturba¬ 
tions  die  out  earlier  than  longer  ones.  This  property  has  been  verified  through  numerical 
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simulations  of  the  type  described  in  the  next  section,  although  we  shall  not  present  evidence 
of  it  in  our  discussion.  Instead,  we  limit  our  attention  to  the  consideration  of  single-mode 
perturbations  only.  This  is  done  in  an  effort  to  underscore  the  differences  in  the  attenuation 
properties  of  the  two  families  of  solutions,  as  well  as  facilitate  comparison  with  our  numerical 
results. 

V.  COMPARISONS  WITH  FAST2D  SIMULATION  RESULTS 

In  this  section,  we  wish  to  test  the  validity  of  Eqs.  (58)  and  (60)  by  comparing  their 
predictions  against  results  from  two-dimensional  computer  simulations.  Several  examples 
are  considered  for  this  purpose  that  illustrate  behavior  from  both  families  of  solutions.  The 
simulations  were  performed  on  a  fixed,  two-dimensional  numerical  grid  using  a  Cartesian  ver¬ 
sion  of  NRL’s  ICF  code  FAST2D  [21],  with  the  thermal  conduction  and  radiation  transport 
modules  turned  off.  Used  in  this  way,  the  FAST2D  code  solves  the  conservation  equations 
of  hydrodynamics  in  Eulerian  form  via  a  flux-corrected  transport  algorithm  (FCT)  [47]. 
All  non-ideal  EOS  data  required  for  this  study  were  derived  from  the  CALEOS  material 
database. 

A  typical  initial  condition  for  our  simulations  appears  in  Fig.  6,  which  shows  a  sinusoidal 
perturbation  superimposed  on  a  two-dimensional  shock  front  moving  into  a  quiescent  homo¬ 
geneous  fluid  medium  on  the  right.  The  size  of  the  computational  grid  was  approximately 
700  x  100  cells,  but  in  Fig.  6,  only  the  first  100  cells  in  the  x-direction  are  shown.  Also  note 
that  this  figure  shows  only  a  perturbed  density  profile,  but  the  pressure  and  x-component 
of  velocity  surfaces  (not  shown)  were  initially  deformed  with  the  same  sinusoidal  ripple. 
The  ^/-component  of  velocity  was  left  unperturbed,  and  initially  set  to  zero  everywhere.  (We 
should  point  out  that  our  choice  of  initial  conditions  here  implies  that  the  Rankine-Hugoniot 
relations  [29]  were  not  strictly  satisfied  in  the  simulations  at  t  —  0;  this  unphysical  situation 
was  quickly  remedied  by  the  FCT  hydro-algorithm  after  a  few  time  steps,  however,  with  no 
noticable  corruption  of  the  numerical  results.)  The  simulations  were  allowed  to  evolve  from 
this  initial  state  for  a  duration  equal  to  at  least  one-and-a-half  ripple-oscillation  periods. 
The  output  from  the  simulations  were  then  post-processed  to  measure  the  evolution  of  the 
ripple  amplitude  as  a  function  of  time. 

A  schematic  of  the  method  used  to  determine  the  temporal  evolution  of  the  shock  ripple 
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amplitude  appears  in  Fig.  7.  The  first  task  was  to  determine  the  position  of  the  deformed 
shock  front  by  computing  a  contour  midway  between  the  unperturbed  upstream  and  down¬ 
stream  density  states.  Since  the  representation  of  a  smoothly-varying  corrugation  to  the 
shock  front  on  a  two-dimensional  Cartesian  grid  is  only  piece-wise  continuous  —  i.e.,  dis¬ 
continuities  exist  between  adjacent  transverse  grid  cells  —  the  initial  contour  in  Fig.  7(a) 
has  a  “stair-step”  appearance,  and  was  Fourier  transformed  and  filtered  to  extract  the  fun¬ 
damental  mode.  The  same  filtering  procedure  was  performed  at  every  subsequent  stage  of 
the  calculation,  as  shown  in  Fig.  7(b)  for  a  time  r  >  0.  The  fundamental  modes  appear  as 
solid  curves  in  the  figures,  along  with  the  positions  of  the  unperturbed  fronts  (solid  vertical 
lines).  The  evolution  of  the  shock  ripple  amplitude  was  then  found  by  computing  —  at  a 
fixed  transverse  location  —  the  distance  between  these  curves  as  a  function  of  time.  Once 
normalized  by  the  initial  shock-ripple  amplitude,  this  distance  (indicated  by  a  double-headed 
arrow  in  each  sub-figure)  gives  an  estimate  for  g(r)/g( 0),  which  can  then  be  compared  to 
theoretical  predictions  based  on  either  Eq.  (58)  or  Eq.  (60). 

Such  a  comparison  is  shown  in  Fig.  8  for  four  different  rippled-shock  systems.  They  are: 
(a)  a  shock  with  M0  =  3  propagating  through  an  ideal  gas  with  7  =  5/3,  where  7  is  the 
ratio  of  specific  heats;  (b)  a  1  Mbar  shock  in  polystyrene;  (c)  a  5  Mbar  shock  in  aluminum; 
and  (d)  a  0.5  Mbar  shock  in  a  cryogenic  mixture  of  deuterium  and  tritium.  (The  Mach 
numbers  in  the  latter  three  cases  were  chosen  so  that  the  internal  energy  of  the  shocked 
material  would  exceed  by  many  times  the  binding  energy  of  the  constituent  atoms,  thus 
warranting  a  hydrodynamic  analysis,  and  justifying  our  use  of  the  term  “fluid  medium”  to 
describe  material  initially  in  a  condensed  state  [29].)  In  Fig.  8,  theoretical  predictions  are 
indicated  by  solid  lines,  and  simulation  results  are  denoted  by  open  circles.  The  relevant 
shock  parameters  for  each  system  are  listed  in  Table  II.  For  the  ideal-gas  shock,  we  see  that 
the  value  of  A  is  positive,  while  the  other  three  examples  belong  to  the  family  of  solutions 
for  which  A  is  negative.  In  all  cases,  the  agreement  between  theoretical  prediction  and 
numerical  simulation  is  quite  good,  which  supports  the  validity  of  Eqs.  (58)  and  (60).  Note 
that  for  the  examples  shown  in  Fig.  8,  the  initial  period  of  oscillation  is  the  time  required 
for  the  shock  to  travel  a  distance  of  about  1  —  2  perturbation  wavelengths  into  the  fluid 
ahead  of  it.  This  period  does  not  remain  constant,  of  course,  but  changes  over  time  and 
asymptotically  approaches  the  value  2nr]/(akD)  —  a  result  that  applies  to  both  families  of 
solutions. 


22 


VI.  SUMMARY  AND  CONCLUSIONS 


In  this  paper,  we  have  generalized  an  earlier  analysis  by  Roberts  [3]  to  derive  explicit 
expressions  governing  the  temporal  evolution  of  perturbations  to  an  isolated  planar  shock 
propagating  through  a  material  with  an  arbitrary  EOS.  It  was  shown  that  under  most 
circumstances,  at  least  two  families  of  stable  solutions  exist.  Membership  in  one  family  or 
the  other  for  a  particular  shock-wave  system  is  determined  by  the  sign  of  the  dimensionless 
quantity  A  [defined  in  Eq.  (55)],  which  is  a  function  of  the  strength  of  the  shock,  and  the 
EOS  properties  of  the  material  through  which  it  propagates.  For  A  >  0,  one  family  of 
solutions  applies  [Eq.  (58)],  while  for  A  <  0,  a  slightly  different  family  [Eq.  (60)]  governs 
the  evolution  of  the  rippled  shock  wave.  Both  families  of  solutions  share  the  same  late-time 
behavior  in  that  the  envelope  of  oscillations  falls  off  asymptotically  as  f”3//2,  but  differ  in  the 
degree  of  damping  that  is  present  initially.  In  general,  solutions  for  which  A  <  0  are  more 
strongly  damped  than  those  with  A  >  0. 

It  is  interesting  to  note  that  the  attenuated  shock-front  oscillations  discussed  in  this  pa¬ 
per  qualitatively  resemble  the  damped  vibrations  of  a  plucked  string  immersed  in  a  viscous 
liquid  [48].  In  the  latter  case,  the  dynamics  are  governed  by  a  differential  equation  sim¬ 
ilar  to  Eq.  (40),  where  the  frictional  term  — ur'(t )  accounts  for  Newtonian  drag  forces  in 
the  surrounding  liquid  that  oppose  the  string’s  motion.  Work  done  against  these  viscous 
forces  by  the  string  drains  kinetic  energy  from  the  vibrating  system  and  converts  it  into 
heat,  resulting  in  oscillations  that  evanesce  over  time.  (We  assume  here  that  the  “radia¬ 
tion  resistance”  of  the  string  [49]  —  which  determines  the  amount  of  energy  converted  into 
sound  —  is  negligible  by  comparison;  since  strings  are  known  to  be  inefficient  radiators  of 
acoustic  energy  [48],  this  approximation  appears  to  be  well  justified.)  As  the  temperature 
of  the  liquid  is  raised,  the  magnitude  of  the  viscous  term  diminishes  [50],  which  lessens  the 
damping  experienced  by  the  string. 

The  analogy  of  a  vibrating  string  in  a  viscous  liquid  provides  some  insight  into  the  dif¬ 
ference  between  the  two  families  of  rippled-shock  solutions  in  Eqs.  (58)  and  (60).  Bearing 
in  mind  the  conventional  microscopic  theory  of  liquids  [51],  this  analogy  suggests  that  the 
strongly-damped  shock-front  oscillations  for  equations  of  state  with  A  <  0  are  likely  a  reflec¬ 
tion  of  appreciable  forces  of  molecular  interaction  in  the  downstream  medium,  particularly 
at  high  densities  and  relatively  low  temperatures  [29,  52],  At  higher  temperatures,  such 
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“viscous”  interactions  become  less  significant,  and  the  behavior  resembles  that  of  a  perfect 
gas  for  which  A  >  0.  These  assertions  are  supported  by  the  fact  that  one  sees  an  even¬ 
tual  (continuous)  transition  from  the  family  of  solutions  with  A  <  0  to  that  with  A  >  0 
for  sufficiently  strong  shocks,  as  shown  in  Fig.  2(c)  for  the  case  of  deuterium-tritium  ice. 
Although  not  indicated  in  Fig.  2,  it  was  noted  during  the  course  of  this  study  that  a  sim¬ 
ilar  sign  change  in  A  occurs  for  shocks  in  polystyrene  and  aluminum  at  approximately  23 
and  50  Mbars,  respectively,  according  to  CALEOS.  [It  should  be  emphasized  here  that  we 
are  not  suggesting  that  the  damping  of  rippled  shocks  is  due  to  the  physical  viscosity  of 
the  downstream  medium,  since  accounting  for  this  fluid  property  lies  beyond  the  Eulerian 
description  adopted  in  Eq.  (12).  Thus,  while  a  vibrating  string  in  a  viscous  liquid  is  a 
suggestive  simplified  model  of  rippled-shock  behavior,  it  should  not  be  taken  too  literally 
in  the  present  context.]  In  the  future,  it  would  be  desirable  to  better  elucidate  the  under¬ 
lying  physical  mechanisms  and  associated  EOS  characteristics  that  are  responsible  for  the 
bifurcated  nature  of  solutions  to  this  class  of  problems. 

The  objective  of  the  present  investigation  was  to  develop  a  better  understanding  of  the 
dynamics  of  rippled  shock  fronts  in  substances  with  non-ideal  equations  of  state.  Stated 
simply,  our  principal  conclusion  is  that  ripple  attenuation  properties  are  EOS  dependent,  and 
can  differ  appreciably  from  those  of  a  perfect  gas,  even  for  moderately  strong  shocks.  This 
result  could  have  important  consequences  for  the  realistic  modeling  of  shock-compressed  ICF- 
fuel  pellets,  since  they  contain  materials  whose  equations  of  state  are  far  from  ideal  (and  often 
poorly  understood).  Because  of  their  potential  for  “seeding”  hydrodynamic  instabilities,  a 
thorough  knowledge  of  how  shock  ripples  evolve  during  the  compression  stage  of  an  ICF 
implosion  is  crucial  for  designing  successful  high-gain  targets.  The  findings  of  this  study 
represents  a  significant  first  step  towards  this  goal.  In  subsequent  investigations,  it  may  be 
possible  to  apply  the  solutions  derived  herein  to  understand  better  the  Richtmyer-Meshkov 
instability  [23]  for  realistic  equations  of  state,  or  to  extend  the  calculation  to  incorporate 
the  influence  of  such  effects  as  convergent  geometries,  and  non-uniform  driving  mechanisms 
( e.g .,  initial  target  roughness  and/or  varying  laser  intensity  in  the  case  of  direct-drive  ICF) 
that  can  launch  perturbed  shocks  [11].  Additionally,  a  more  complete  understanding  of  the 
dynamics  of  isolated  rippled  shocks  may  be  useful  in  the  study  of  certain  type  II  supernova 
phenomena  [53] ,  and  forms  the  basis  for  developing  a  new  analytical  benchmark  to  validate 
the  performance  of  ICF  and  “high  energy-density  physics”  codes  [17-20]. 
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APPENDIX:  LAPLACE  TRANSFORMS  AND  MATHEMATICAL  IDENTITIES 


We  cite  in  this  appendix  particular  Laplace-transform  relations  and  miscellaneous  math¬ 
ematical  identities  that  were  useful  in  our  analysis  of  the  isolated  rippled-shock  problem. 
Following  convention,  we  let  the  symbol  £  represent  the  Laplace  transform  of  a  function 
<3>(t),  defined  as 

POO 

£{$(r)}=  /  e-ST$(r)dr  =  $L(S),  (Al) 

Jo 

where  the  subscript  L  is  a  practical  shorthand  notation,  and  s  is  the  Laplace  transform 
variable.  The  inverse  transform  operator  is  denoted  by  /A1,  such  that  £_1  ($l(s)}  =  <3>(t)  . 
In  terms  of  these  definitions,  the  following  results  can  be  established  [54] : 

£{1}  =  1/s  ,  (A2) 


£{J0(ar)}  =  l/Vs2  +  a2  , 


(A3) 


£ 


-i 


c 


-1 


£  { J\(a  r)}  =  a  (  1 


£  <(  /  Jo(w)  dw 
>o 


s  +  Vs2  +  a2 

1 

(s-n)2  +  n2 


Vs2  +  a2  J  ’ 
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iVs2  + 


ar 


Ji(ar) 


ar 


= 


sin  Qt 

n  ’ 


(A4) 

(A5) 

(A6) 

(A7) 
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-i 


——J——  j  =  (cos  ^  sin  Qr 


(AS) 


Note  that  the  constants  a,  Q,  and  //  appearing  in  the  equations  above  are  real  quantities.  An 
important  consequence  of  Laplace-transform  theory  is  the  convolution  theorem  [41],  which 
asserts  that  any  two  continuous  and  sufficiently  “well-behaved”  functions  <h(r)  and  T(r) 
obey  the  relation 

£-1{$l(s)Tl(s)}  =  [  $(T-z)*(z)dz.  (A9) 

Jo 

Stated  differently,  the  convolution  theorem  says  that  the  product  <Ll(s)'Fl(s)  is  the  Laplace 
transform  of  the  function  defined  by  the  right  side  of  Eq.  (A9)  —  a  result  that  plays  a  central 


role  in  the  derivation  of  the  solutions  presented  in  Sec.  IV.  Other  important  mathematical 
identities  for  this  study  were  [55] : 


(  VA1 2 3 4 5 6 7 8  -  a2  -  A 


Ji(az) 


cr 


cos  At  if  A  >  a, 


az 


sin  A (r  —  z)dz  =  < 
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- - -  sm  At - -  cos  at  if  A  <  a, 
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where  A,  cr,  and  a  are  positive  real  parameters. 
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Ideal  gas  (7  =  5/3,  Mq  =  3) 


CH  (1  Mbar) 


A1  (5  Mbar ) 


DT  (0.5  Mbar ) 


PO 

— 

1.07  g/cm 3 

2.71  g/cm 3 

0.25  g/cm3 

pi 

3po 

2.350  g/cm3 

6.257  g/cm 3 

0.7903  g/cm3 

T0 

— 

300  K 

300  K 

19  K 

T, 

11  T0 

5.290  x  10 3  A 

1.879  x  10 4  A 

8.969  x  103  I< 

Po 

__ t 

1.779  x  10-2  Mbar 

1.813  x  10-2  Mbar 

8.917  x  10“4  Mbar 

Pi 

11  Po 

1  Mbar 

5  Mbar 

0.5  Mbar 

D 

V15Po/Po 

1.298  x  106  cm/s 

1.801  x  106  cm/ s 

1.709  x  106  cm/s 

c 

\/55po/{9po ) 

1.236  x  106  cm/s 

1.523  x  106  cm/s 

1.207  x  106  cm/s 

rj 

3 

2.197 

2.309 

3.161 

h 

-0.1111 

-0.1647 

-0.2071 

-8.448  x  10”2 

Mo 

3 

5.868 

3.355 

11.72 

Mi 

0.5222 

0.4782 

0.5120 

0.4478 

a2 

3.577 

2.006 

1.983 

3.413 

0? 

2.667 

3.375 

2.812 

3.988 

a\ 

7.289 

7.343 

6.449 

12.24 

P 

1.064 

1.218 

1.179 

1.211 

r 

2.553 

1.919 

1.788 

3.232 

a 

2.961 

1.433 

1.527 

2.103 

b 

-1.316 

-0.8412  i 

-0.7303  i 

-0.5551* 

(< a  +  b )2 

2.706 

1.346  —  2.411  i 

1.799-  2.230  i 

4.115  —  2.335  * 

(a  —  6) 2 

18.29 

1.346  +  2.411  i 

1.799  +  2.230  * 

4.115  +  2.335  * 

a 

—1.316  i 

0.8412 

0.7303 

0.5551 

A 

4.214 

-5.419 

-3.015 

-4.747 

T  Assuming  po  and  To  are  independent,  we  have  po  =  1Z  poTo/m,  where  1Z  =  8.317  x  10'  erg/deg  • 
mole  is  the  universal  gas  constant,  and  m  is  the  molecular  weight  of  the  gas. 


TABLE  II:  Parameters  for  the  isolated  rippled-shock  examples  considered  in  this  paper.  The 
labels  CH ,  Al,  and  DT  stand  for  polystyrene,  aluminum,  and  deuterium-tritium,  respectively. 
The  entries  in  the  last  three  columns  were  computed  using  the  CALEOS  database. 
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FIG.  1:  Qualitative  difference  between  the  two  families  of  solutions  for  an  isolated,  rippled  shock 
wave.  Although  the  amplitudes  of  shock  ripples  decay  asymptotically  in  both  cases  as  f-3/2,  the 
degree  of  initial  damping  in  one  family  is  typically  smaller  than  in  the  other.  The  lesser-damped 
example  shown  in  (a)  corresponds  to  a  positive  value  of  the  parameter  A  [defined  in  Eq.  (55)], 
whereas  in  (b),  A  is  negative. 
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FIG.  2:  Examples  of  Hugoniot  curves  for  (a)  polystyrene,  (b)  aluminum,  and  (c)  deuterium-tritium 
ice.  Solid  and  dashed  curves  show  results  derived  from  the  CALEOS  and  SESAME  EOS  libraries, 
respectively.  The  symbols  in  (b)  are  experimental  data  points  taken  from  Ref.  [28] .  According  to 
CALEOS,  the  sign  of  A  in  (c)  changes  from  negative  to  positive  near  p±  =  1.7  Mbar,  indicating  a 
transition  to  a  lesser  damped  solution  above  that  value  (shaded  region).  For  the  range  of  density 
values  in  (a)  and  (b),  the  sign  of  A  is  strictly  negative. 
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FIG.  3:  Geometry  of  an  isolated,  rippled  shock  front  in  the  frame  of  reference  in  which  the  down¬ 
stream  fluid  is  at  rest.  The  unperturbed  planar  shock  moves  in  the  negative  x-direction  with 
velocity  —(D  —  U )  x,  while  the  upstream  fluid  flows  at  velocity  U  x.  Unit  vectors  normal  and 
tangential  to  the  rippled  front  are  denoted  by  N  and  T,  respectively. 
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FIG.  4:  Space-time  diagram  of  shock  and  sound  wave  that  define  the  solution  region  G  at  the  time 
ra ■  This  region  is  bounded  by  the  countour  C,  which  is  comprised  of  three  line  segments:  Ci, 
Ci i,  and  Cm.  The  first  extends  along  the  £-axis  from  0  to  ra(M^ 1  —  1),  the  second  lies  on  the 
characteristic  £  +  ra  =  (ra  —  r) /Mi,  and  the  third  coincides  with  the  path  of  the  shock  wave  given 
by  £  +  t  =  0. 
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A 


FIG.  5:  Plot  of  A  =  a4  —  4f5Ta2  +  4T2  as  a  function  of  a2 .  The  roots  of  A  =  0  are  denoted  by 
a 2_  and  as  indicated  in  the  figure.  If  A  <  0,  the  value  of  a 2  lies  between  these  two  limits.  For 
physical  equations  of  state  with  A  >  0,  the  value  of  a2  is  apparently  restricted  to  be  less  than  a2  . 
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FIG.  6:  An  example  of  a  perturbed  density  surface  at  the  start  of  a  FAST2D  simulation.  The 
ratio  of  ripple  amplitude  to  wavelength  in  this  study  is  5%.  Periodic  boundaries  are  assumed  in 
the  transverse  direction,  and  inflow  and  outflow  conditions  are  imposed  at  the  left  and  right  ends, 
respectively,  of  the  computational  domain. 
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FIG.  7:  Schematic  of  the  method  used  for  determining  the  shock  ripple  amplitude  in  the  FAST2D 
simulations  (a)  initially,  and  (b)  at  a  later  time  r  >  0. 


FIG.  8:  Comparison  of  theoretical  predictions  (solid  lines)  and  FAST2D  simulation  results  (open 
circles)  for  the  normalized  ripple  amplitude  of  a  perturbed  shock  wave  propagating  through  (a) 
an  ideal  gas,  (b)  polystyrene,  (c)  aluminum,  and  (d)  deuterium-tritium  ice.  In  (a),  the  ratio  of 
specific  heats  is  7  =  5/3  and  the  unperturbed  Mach  number  is  Mq  =  3.  In  (b),  (c),  and  (d),  the 
unperturbed  shock  strengths  are  1,  5,  and  0.5  Mbar,  respectively.  The  example  in  (a)  belongs  to 
the  family  of  solutions  for  which  A  >  0;  those  in  (b),  (c),  and  (d)  correspond  to  A  <  0. 


